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ABSTRACT 


An  analytic  semi-empirical  approximation  to  the  extinction  efficiency,  Qext,  for 
randomly  oriented  spheroids,  based  on  an  extension  of  the  anomalous  diffraction  formula,  is 
given  and  compared  with  the  extended  boundary  condition  or  T-matrix  method.  This  will 
allow  for  better  and  more  general  modelling  of  obscurants.  Using  this  formula,  Qext  can  be 
evaluated  over  104  times  faster  than  with  previous  methods.  This  approximation  has  been 
verified  for  complex  refractive  indices  m  =  n  —  ik,  where  1  <  n  <  oo  and  0  <  k  <  oo  and 
aspect  ratios  from  0.2  to  5.  We  believe  that  the  approximation  is  uniformly  valid  over  all 
size  parameters  and  aspect  ratios.  It  has  the  correct  Rayleigh,  refractive  index  and  large 
particle  asymptotic  behaviours.  The  accuracy  and  limitations  of  this  formula  are  extensively 
discussed. 


RESUME 

Nous  presentons  une  approximation  analytique  semi-empirique  de  l’efficacite 
d’extinction  Qext,  pour  des  particules  spheroi’dales  orientees  aleatoirement,  basee  sur  une 
extension  de  la  theorie  de  la  diffraction  anomalique.  Cette  approximation  est  com’  ee 
au  calcul  base  sur  la  methode  de  conditions  frontieres  etendues  ou  matrices  de  tran  a. 
Cette  approximation  permet  de  mieux  modeliser  les  obscurants.  Grace  a  cette  formule,  ^ert 
peut  etre  evalue  104  fois  plus  vite  que  par  les  method es  anterieures.  Nous  avons  verifi4  cette 
formule  pour  des  indices  de  refraction  complexes  m  =  n  —  ik,  oul<n<ooetO<Jb<oo 
et  pour  des  aspects  de  0.2  a  5.  Nous  croyons  que  ^approximation  est  uniformement  valide 
pour  toutes  les  grosseurs  et  tous  les  aspects.  Elle  se  comporte  correctement  a  la  limite  de 
Rayleigh,  a  la  limite  de  l’indice  de  refraction  infini  et  a  la  limite  geometrique.  L’exactitude 
de  cette  formule  ainsi  que  ses  limites  sont  discutees  en  detail. 
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The  present  and  future  electro-optic  capabilities  of  fire  control  systems  require  full 
spectrum  obscurants  as  an  effective  countermeasure.  It  is  well  known  that  spherical  particles 
can  obscure  efficiently  only  in  the  visible  region  of  the  spectrum.  In  the  near  and  far  infrared 
as  well  as  in  the  microwave  region  of  the  spectrum  the  use  of  non-spherical  particles  is 
mandatory.  Imaging  through  ocean  waters  and  dusty  environments  both  require  knowledge 
of  the  extinction  by  non-spherical  particles.  Communications,  remote  sensing,  surveillance 
and  radar  propagation  frequently  encounter  scenarios  involving  non-spherical  particles.  The 
estimation  of  electromagnetic  extinction  through  dust  or  other  non-spherical  atmospheric 
aerosols  and  hydrosols  is  an  essential  first  step  in  the  evaluation  of  the  performance  of  all 
electro-optic  systems. 

To  date,  the  theoretical  exploration  of  the  effect  of  these  types  of  particles  on  the 
performance  of  obscurants  has  been  either  extremely  restrictive  or  prohibitively  expensive. 
This  work  partially  overcomes  these  severe  restrictions  and  will  aid  in  the  design  and  perfor¬ 
mance  analysis  of  potential  new  obscurants.  For  instance,  a  factor  greater  than  ten  thousand 
in  computational  speed  has  been  achieved  with  minimal  loss  in  accuracy.  Furthermore,  the 
present  approach  is  much  more  flexible  in  the  sense  that  extremely  elongated  or  flattened 
particles  can  now  be  modelled.  Elongated  particles  are  thought  to  be  good  millimeter  wave 
obscurants  and  flattened  or  flake-like  particles  are  known  to  be  good  infrared  obscurants. 

The  immediate  objective  of  this  work  is  to  significantly  reduce  the  computational 
burden  in  calculating  the  extinction  from  non-spherical  particles.  This  allows  for  the  explo¬ 
ration  of  the  effects  of  non-spherical  particles. 

The  longer  term  goal  of  this  work  is  to  alleviate  the  remaining  constraints  in  the 
theoretical  consideration  of  non-spherical  aerosols  and  obscurants  of  either  natural  or  ar¬ 
tificial  origin.  This  will  aid  not  only  in  finding  better  obscurants  but  also  in  the  possible 
identification  and  remote  classification  of  such  aerosols. 
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1.0  INTRODUCTION 

The  immediate  objective  of  this  work  is  to  significantly  reduce  the  computational 
burden  in  calculating  the  extinction  from  non-spherical  particles.  This  allows  for  the  ex¬ 
ploration  of  the  effects  of  non-spherical  particles  on  the  performance  of  obscurants  and 
electro-optical  systems.  The  longer  term  goal  of  this  work  is  to  alleviate  the  remaining  con¬ 
straints  in  the  theoretical  consideration  of  non-spherical  aerosols  and  obscurants  of  either 
natural  or  artificial  origin.  This  will  aid  not  only  in  finding  better  obscurants  but  also  in 
the  possible  identification  and  remote  classification  of  such  aerosols. 

We  have  previously  presented  (Ref.  1)  a  numerical  approximation  to  Qext  for  ran¬ 
domly  oriented  spheroids.  This  work  was  applicable  to  particles  with  1.01  <  n  <  2  and 
0  <  k  <  1  for  arbitrary  sizes  and  aspect  ratios.  The  required  angular  integration  was 
carried  out  by  a  64-point  Gaussian  quadrature.  Since  many  materials  have  optical  prop¬ 
erties  beyond  the  above  limits,  we  have,  in  this  work,  extended  the  refractive  index  range 
to  1  <  n  <  oo  and  0  <  k  <  oo.  Furthermore,  large  optical  sizes  produce  high  frequency 
oscillations  in  the  kernel  of  the  angular  integral.  These  integrals  are  very  difficult  and 
time-consuming  to  estimate  numerically.  We  have  replaced  this  numerical  integral  with  an 
approximate  analytic  expression  that  overcomes  this  difficulty. 

The  basic  approach  is  to  orthogonalize  as  much  as  possible  the  scattering  physics 
into  well-defined  regimes.  For  small  physical  and  optical  sizes,  the  electrostatics  (Rayleigh) 
approximation  is  used.  For  larger  and  very  large  optical  sizes  we  still  use  the  electrostatics 
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approximation  but  with  the  optica]  constants  transformed  (Ref.  2)  to  include  the  effects  of 
the  magnetic  dipole.  In  the  large  physical  si2e  regime  we  split  the  physics  into  a  diffraction 
(anomalous  diffraction)  component  and  what  can  be  loosely  described  as  edge  effects.  The 
diffraction  component  is  modelled  by  the  anomalous  diffraction  approximation  as  developed 
by  Van  de  Hulst  (Ref.  3).  The  edge  effect  (Fock  theory)  component  is  modelled  by  extending 
a  technique  introduced  by  Jones  (Refs.  1  and  4).  In  this  report,  this  component  is  further 
generalized  to  have  proper  behaviour  at  small  optical  sizes  and  for  large  indices.  The  gap 
between  the  large  and  small  particle  regimes  is  bridged  by  a  binomial  form  (Ref.  1)  similar 
to  the  generalized  mean  (Ref.  5,  p.  10). 

The  report  is  organized  as  follows:  Chapter  2.0  develops  the  extinction  formula. 
This  includes  Sections  2.1,  small  particle  scattering,  2.2  anomalous  diffraction,  2.3  the  edge 
effects  and  2.4  the  binomial  bridging  function.  Chapter  3  and  the  Appendix  contain 
extensive  comparisons  of  this  approximation  to  the  T-matrix  method  (Ref.  6).  The  final 
chapter,  Chapter  4.,  gives  the  conclusions  and  remaining  limitations. 

This  work  was  performed  at  DREV  between  October  1990  and  March  1992  under 
PSC  32A,  EO/IR  Protection  of  Land  Vehicles. 

2.0  DEVELOPMENT  OF  THE  EXTINCTION  FORMULA 

The  physics  and  the  structure  of  the  extinction  formula  for  randomly  oriented 
spheroids  will  be  developed  in  this  chapter. 
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2.1  Small  Particle  Term 

For  a  physically  and  optically  small  particle,  all  applied  electric  field  gradients 
disappear  and  the  particle  begins  to  respond  to  a  homogeneous  field.  The  electrostatics 
approximation  then  holds,  giving  rise  to  the  Rayleigh  scattering  formula  for  randomly  ori¬ 
ented  spheroids  (Ref.  7).  If  the  particle  is  still  geometrically  small  but  optically  large,  the 
magnetic  dipole  field  becomes  significant.  We  have  found  (Ref.  2)  that  one  can  use  the 
electrostatics  approximation  but  with  the  optical  constants  transformed.  This  transform 
correctly  describes  the  full  field.  It  is  exact  for  spheres  and  approximate  for  randomly  ori¬ 
ented  spheroids.  The  following  formulaefor  the  extinction  efficiency  from  small  particles, 
Qtmaii  axe  derived  in  Ref.  2: 


Q  small  —  Q  sea  "f  Qabs 


with 


where 


=  y^{lm|!  +  IS|2+2(lfe|!  +  l5|1)}  and 

<3^.  =  +  + 2  (%  +  >£))) 


Vi  = 


3(X*  +  ?IPt)  ’  V'  *(Li  +  wiTTl)  ’ 


[1] 

[2] 

[3] 


The  optical  constants  transform  is 
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where  *l>\(z)  is  the  first  Ricatti- Bessel  function  and 


Z\  =  y/cjib(l  +  x(l  -  1/r2)) 

prolates 

= 

oblates 

[5] 

z2  =  y/eJib(l  +  x1/3^-^/y/r)) 

prolates 

=  y/ejibi^'3 

oblates 

[6] 

where  x  ~  >/3/10  and  £  «  {(91|m  -  1|  -  64)/155  |  0  <  C  < 

!}• 

The  normalization  factor  is 

_  r2  ,  /  s/t1  —  1  \ 

A  =  1+  r-a, — "  s‘n  -  for  prelates, 

v^nri  \  r  J 

r2  / 1  4-  Vl  —  r2\ 

=  1  +  —  In  ( - )  for  oblates. 

Vl  -  r2  \  r  J 
The  form  factors  are  defined,  for  prelates  as 

Li  =  — -2— +  )},  L2  =  ~2-,  5  =1-^ 


[7] 


[8] 


and  for  oblates, 


Lx  = 


L2  =  ^I,  /-l-l. 


l±/2fi  tan-1/ 

P 


[9] 


The  above  set  of  equations  is  identical  to  that  found  in  Ref.  2  apart  from  the 
exponent  C  in  [4].  This  exponent  is  introduced  for  the  following  reason.  The  exact  solution 
of  the  scattering  problem  can  be  expressed  in  terms  of  scattering  coefficients  (for  spheres, 
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these  are  the  Mie  coefficients).  For  small  physical  sizes  the  first  scattering  coefficient  is 
dominant  and  can  be  expanded  as  an  infinite  power  series  in  the  size  parameter  x.  The 
first  term  of  this  series  is  dominant  only  when  |m  -  1|  is  small.  If  we  transform  the  optical 
constants  then  the  first  scattering  coefficient  can  still  be  expressed  as  an  infinite  series  in  x, 
but  this  time  the  first  term  dominates  only  when  |m—  1|  is  large.  For  intermediate  values  of 
|m  —  1|,  more  terms  of  the  power  expansion  of  the  scattering  coefficients  are  required.  Since 
the  resulting  expression  would  be  cumbersome  we  have  decided  to  modify  the  first  term 
of  the  transformed  series  in  such  a  way  as  to  model  empirically  the  scattering  behaviour 
from  small  to  large  |m  —  1|.  A  simple  and  robust  way  to  achieve  this  is  to  gradually  turn 
on  the  transform  by  using  the  exponent  (.  Thus  when  C  =  0  the  Rayleigh  expression 
results  and  when  C  =  1  the  fully  transformed  expression  results.  Note  that  C  =  0  when 
|m  —  1|  <  64/91  at  0.7  and  C  =  1  when  ]m  -  1|  >  219/91  at  2.4. 

2.2  Large  Particle  Term 


In  this  section  we  discuss  the  extinction  efficiency  of  large  particles  Q iarge.  We 
separate  the  physics  into  two  parts,  one  which  corresponds  to  the  anomalous  diffraction, 
Qad ,  and  the  other,  which  can  be  considered  due  to  edges  effects,  Qed- 


UNCLASSIFIED 

6 


2.2.1  Extended  Anomalous  Diffraction 


The  anomalous  diffraction  formula  is  derived  (Refs.  1  and  8)  by  assuming  that  the 
incident  plane  wave  is  not  significantly  skewed  in  passing  through  the  scattering  object  and 
that,  to  first  order,  the  effect  of  the  scatterer  is  to  locally  delay  the  phase  of  the  wave  and 
attenuate  its  amplitude  (Ref.  3).  The  strict  limit  of  validity  of  the  formula  is  therefore 
the  region  where  (n  -  1)  <<  1.  The  scattering  object  is  in  effect  treated  as  an  irregular 
disc  normal  to  the  incident  wave  possessing  a  spatially  dependent  phase  and  amplitude. 
The  Fraunhofer  pattern  at  infinity  is  then  derived  and  Qext  evaluated  from  the  standard 
relations.  For  a  spheroidal  scatterer,  this  procedure  leads  immediately  to  the  following 
formula  : 

Q.,  =  R«{2  +  4£l  +  4£9^}  [10] 


where  u  is  given  by 


ui  =  »Ai/>, 


A i/>  =  2  (m  -  1)— 
P 


[11] 


and 


p  =  V cos2  6  +  r2  sin2  0,  a  =  2xa/A,  b  =  2x/?/A,  m  =  n  -  ik.  [12] 


Where  r  =  a/b  is  the  aspect  ratio  (for  prolates  r  >  1  and  for  oblates  r  <  1),  a  is  the  length 
of  the  semi-axis  of  rotation,  /?  is  the  other  axis  of  the  spheroid,  0  is  the  angle  between  the 
incident  radiation  and  the  a  or  a  axis,  and  X  is  the  wavelength  of  the  scattered  radiation. 
Hence  a  and  6  are  the  two  size  parameters  associated  with  the  spheroid,  and  p  can  be 
considered  a  projection  operator  of  the  penumbral  ellipse  (the  ellipse  defined  by  the  shadow 
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line  on  the  surface  of  the  spheroid)  onto  the  plane  perpendicular  to  the  direction  of  the 
incoming  radiation.  When  compared  with  exact  results  this  formulation  is  satisfactory 
for  all  oblate  spheroids.  For  prolates,  however,  significant  phase  differences  appear  due  to 
deviation  of  the  central  ray. 


To  account  for  the  refraction  of  the  central  ray,  we  have  previously  modified  the 
anomalous  diffraction  approximation  (Ref.  1).  The  deviation  of  the  central  ray  is  now 
taken  into  account  when  computing  its  phase  difference. 


The  extended  anomalous  diffraction  A ij>  is  found  to  be 


w  =  iAip 


cos(<f>) 

sin(<f>) 

A 


.,  f2rr _ P2 cos (<ft)  +  s sin(<ft) _ ll.  _ 

[pip2  cos2(<£)  +  q 2  sin2(<£)  +  2s  cos(<£)  sin (<j>)  if™ 

s2  +  p2  A 


m(p4  +  s2) 

-  p2) 

m(p4  +  s2) 


[m2(p4  +  52)-s2]1/2 


q  =  [r2cos2(«)  +  sin2(«)]1/2.  [13] 

In  the  limiting  cases  of  r  — ►  oo  A0  becomes: 


A^>  =  26{(m2  -  cos 2  Off2  -  sin  0). 


[14] 


For  random  orientations  the  angular  averaging  is  carried  out  as  follows: 

q  fo /2QadP  sin  0d0 
fo'2  p  sin  0d0 


[15] 
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The  integration  in  the  numerator  of  [15]  can  be  readily  computed  numerically  if 
the  kernel  is  not  too  oscillatory.  However,  [15]  can  be  analytically  approximated  by  the 
technique  described  in  Ref.  9.  This  eliminates  numerical  difficulties  and  leads  to  a  more 
efficient  algorithm.  Following  the  procedure  of  Ref.  9,  [15]  becomes: 

Qad  =  2  +  4(h-h)/m  [16] 


with 


h  =  A  ]e~c 


{<1+e> 


C2 


C2 


[17] 


and 


where 


/2  = 


m  a  i 

C2  r  C2M 0) 


1  x  ,  2A  Bu(x/2)  1 

w(jt/2  )'r  C3  [(B  +  j{0))u(0). 


[18] 


A 

C 


7  B  [BA  j(0)], 
w(0)-7  B, 


_ u(0)  -  u>(7t/4) _ 

7  +  [w(ir/4)  -w(ir/2)]/j(ir/4)’ 
u/(tt/2)  —  w(0) 

i(0)  ’ 


[19] 


and 


m 


1 

2 


cos(#)\/l  -  g2cos2(fl)  +  sin~I(^cos(^))  f 
co&(0)y/l  +  P  cos(0)2  +  M/cosig)iVU^coS(^) 


for  prelates, 


for  oblates.  [20] 


Also, 


F  fn  =  En{u(0)  ~  c)  -  W*/2)  -  g) 

M  W  ~  (u>(0)  -  C)"-1  Htt/2)  -  C)"'1  ’ 


[21] 
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where  En  is  the  nth  order  exponential  integral  (Ref.  5).  And  finally  for  prolates  (considering 
deviated  rays), 


«(0) 

w(ir/4) 


2b  i 
2b  i 


\/m2  —  1  (m  -  1) 

(1  -  a)(m  -  1)  +  aVm2  -  1/r 


[l72(\/m2  -  1/2  -  l/\/2)  +  . 


w(jt/2) 

a 


2 bi(m  -  1), 


1  -  rc 


(1-  0 

(1  —  rc)’ 


rc  = 


/ 

m  -  1 

V 

m  +  1 

for  oblates  (undeviated  rays), 


[22] 


u;(0)  =  26ri(m-l), 

u{ir/A)  ~  2bri(m~ 

w(ir/ 2)  =  2bi(m  -  1).  [23] 


As  a  consequence  of  the  analytic  integration  technique,  [13]  is  approximated  by  [22] 
at  three  values  of  6  whereas  [11]  exactly  reduces  to  [23]  at  the  same  three  values.  Note  that 
[16]  can  be  shown  to  reduce  to  the  anomalous  diffraction  formula  of  Van  de  Hulst  for  the 
sphere. 


2.2.2  Edge  Effects 


For  a  particle  whose  typical  size  is  much  larger  than  the  wavelength,  the  edge 
cannot  be  treated  as  sharp  and  the  effect  of  the  curvature  of  the  object  must  be  included. 
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Jones  (Ref.  4)  has  shown  how  to  estimate  these  edge  effects  for  three  dimensional  convex 
bodies.  In  Ref.  1  we  showed  that 


It  can  be  shown  for  convex  bodies,  randomly  oriented  or  illuminated  by  a  randomly  polarized 
beam,  that  ct  is  a  universal  function  of  m.  As  |m  -  1|  — ►  0,  ce  — ►  cq  =  0.996130  and  as 
|m  -  1|  — >  oo,  ce  -*  Cqo  =  0.0659708  (Refs.  1,  11  and  12).  This  universal  function  is 
not  known  and  therefore  must  be  modelled.  Since  ce  is  independent  of  shape,  this  can  be 
accomplished  by  studying  the  sphere  alone.  An  added  complication  to  D  arises,  however, 
when  we  consider  spheroids  with  small  phases,  i.e.  A^>  <<  1.  This  occurs  since  what 
we  have  been  calling  an  edge  effect  is  in  fact  the  field  distortion  around  the  boundaries 
of  the  particle,  and  hence  its  behaviour  for  small  A  ip  is  quite  different  than  for  large  A^>. 
As  before,  we  have  modelled  this  effect  in  our  expression  for  D  by  using  the  sphere.  Our 
empirical  model  of  the  above  two  behaviours  is: 


[6W/3  -I-  l/|m  -  lp]V0 

_ 4/25 _  C0  —  Cqq 

|m  — 1|3  +  8/125  In  |m  -  lj  ’ 


|mj  <  1.95  or  |m|  >  2.05 


1.95  <  |m|  <  2.05 


[26] 
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In  the  above,  the  special  case  for  f3  with  m  as  2  is  to  avoid  the  obvious  singularity  in  the 
maun  expression.  At  this  point  we  require  only  a  relatively  large  value  for  /?. 


The  large  particle  limit  for  the  spheroid  becomes 


Qext  *  2  +  Qedge • 


[27] 


We  now  wish  to  produce  a  term  T  which,  when  it  multiplies  [10],  gives  the  same  limit  as 
[27]  without  diverging  as  the  size  parameter  goes  to  zero.  We  have  found  that  (Ref.  1)  an 
adequate  expression  for  our  purposes  is: 


T  =  2-  e-Q‘d»‘/2 


[28] 


For  random  orientations  the  angular  averaging  is  carried  as  in  [15]  but  with  T  replacing 

Qad ‘ 

T  =  j^Tp.iuide  [29J 


f^2  p  sin  Odd 

Using  the  same  integration  technique  to  obtain  am  analytic  approximation  to  [29] 
we  obtain: 

T  =  2  -  e~sc'  A'F^/j'(0)  [30] 


where 


a'  =  yu'[*'+;'(o)],  s'  =  7TR 

c  = 


1  -  v(r/4) 


x/4)  -  v(*/2)]/i'(ir/4)’ 


1  _  y  b'  V  =  K*/2)  ~  1  $  -  dt 2/3 

1  7  °  ’  7  -  j>( 0)  ’  0  -  UT  • 


and 


,  _  W)  _  E2(6v(x/2)) 
6  6v(ir/2) 


[31] 


[32] 
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Here  £3,  is  the  2nd  order  exponential  integral.  For  prelates, 


m  =  ' m. 

«(»/4)  =  1f,[-!,i;l;(r!-l)/(-'!  +  l)|(2/(rJ  +  l)|3 

»(t/2)  =  [33] 

and  for  oblates, 

m  =  - as ), 

v(jt/4)  =  2£il-|,^;l;l-r2/2]/r 

«(jr/2)  =  2^i[-|,^;l;l-r2]/r2  [34] 

2.2.3  Total  Contribution  to  QiaT3e 


The  total  contribution  to  Q  large  for  a  given  orientation  is  (Ref.  1): 

Q  large  = 


[35] 


.or  random  orientation  the  angle  averaging  would  give: 

fo /2QadTp  sin  6d6 
p sin  Odd 

Due  to  their  complementary  nature,  Qa 4  as  2  =  constant  for  large  |u;|  while  T  as  constant 
for  small  to  medium  |ur|.  Therefore,  to  a  good  approximation,  we  can  separate  the  kernel 
of  the  above  integral  to  obtain: 


<2< 


large 


^2  large  =  ^2  ad?  • 


[37] 
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An  extensive  comparison  between  [37]  and  the  numerical  computation  of  <7/orfle  from  [36] 
provides  definitive  evidence  for  the  above  argument.  Some  examples  of  this  are  shown  in 
the  next  chapter. 


2.3  Bridging  Function 


^ small  *s  a  g°°d  approximation  to  £?eit  when  the  semi-major  axis  is  less  than  1. 
For  larger  values  of  the  semi-major  axis  it  overestimates  Qexi-  Similarly,  QtaTge  is  a  good 
approximation  to  £?exl  when  the  semi-major  axis  is  greater  than  2-5  depending  on  |m—  1|. 
For  smaller  values  of  the  semi-major  axis  it  overestimates  To  obtain  $ext  in  the 

transition  region  from  ($iarge  and  Qtmalh  we  need  a  bridging  function  that  smoothly  goes 
between  the  two.  The  form  must  have  Ismail  88  *be  brst  term  in  its  series  expansion  and 
asymptotically  go  to  ~Qlarge. 


A  quite  general  form  that  can  do  this  is  the  confluent  hypergeometric  function  or 
Kummer  function  which  has  the  general  form  iFi[a-,b,cz''],  where  a,b,c,v  are  arbitrary 
parameters  and  z  is  tne  variable  (Ref.  5).  From  the  basic  properties  of  this  function 


lim  r(6  -  a) 

*—•  oo  r(6) 


cazav 


c*za\F1[a]b,-cz1'] 


£ (b  -  a) 

m 


l  +  o(l±a-6)  + 
czv 


{*-¥ 


2u 


ba(a  -|- 1  )z 

2!  6(6  +  1) 


-} 


[38] 


With  the  small  and  large  particle  limits  considered  the  bridging  function  B  becomes: 


B  =  Qsmcirf'i[l/v;ii.~{czy\  with 
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c  =  r(6)/r(6  - 1/|/),  and  z  =  ’Qmall/^largt  [39] 

which  has  the  correct  limits,  i.e. 


limB  -  ~Q,maU 

z—+0 

}^oB  ^ large  [40] 


We  now  need  to  determine  v  and  b.  We  should  expect  these  parameters  to  be  insensitive  to 
aspect  ratio  since  most  of  the  shape  effects  are  already  accounted  for  by  and  (?/orje. 

We  have  confirmed  this  supposition  by  numerical  evidence.  Assuming  these  parameters  are 
shape  independent,  we  can,  for  each  of  an  array  of  m  values,  find  the  ‘best’  set  of  values 
of  v  and  6  by  considering  the  sphere  only.  We  did  this  by  using  a  nonlinear  fitting  routine. 
The  results  of  this  fit,  for  the  vast  majority  of  cases  made  b  very  large  (>  10).  When  b 
becomes  very  large,  B  goes  confluent  (see  Ref.  13  for  more  details).  This  means  that  the 
function  simplifies  dramatically  to  a  binomial  function,  thus: 


lim  B 

b— *oo 


^ email 

[1  +  z^/u 


This  can  be  rearranged  in  a  form  similar  to  a  generalized  mean  as: 


[41] 


large 


[42] 


We  must  now  model  the  values  for  7.  We  have  modelled  this  parameter  in  previous  papers 
(Refs.  1  and  14).  This  previous  modelling  is  not  useful  to  the  current  approach  since  ($^011 
and  iarge  have  changed,  the  formula  for  ^ext  has  been  extended  to  all  n  >  1  and  angle 
averaging  has  been  carried  out. 
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In  [42]  if  one  of  either  (5,maii  or  (Jiorfle  is  much  smaller  than  the  other  then,  to  a 
good  approximation,  is  this  value  provided  that  7  >  0.  Larger  values  of  7  will  drive 
to  the  smallest  value  between  and  C J/orje  more  than  smaller  values.  When  the 

physical  size  of  the  particle  is  small,  Qmau  is  equal  to  or  less  than  Qiarge  and  hence  7  should 
be  large.  Conversely,  when  the  physical  size  of  the  particle  is  large,  Qmaii  is  much  larger 
than  Q iarge  and  we  then  merely  require  that  7  be  positive.  The  real  difficulty  in  modelling 
7  is  thus  in  the  intermediate  ranges  of  physical  sizes.  Most  of  the  behaviour  of  7  in  the  large 
particle  regime  can  be  determined  by  considering  first  the  real  axis,  in  the  refractive  index 
plane,  and  then  by  varying  the  imaginary  component.  Considering  first  the  real  axis,  the 
behaviour  can  be  described  by  the  sum  of  two  asymptotic  terms,  one  as  n  -  1  — ►  0  and  the 
other  as  n  — >  00.  The  behaviour  of  7  as  A:  increases  i6  well  modelled  by  a  single  term.  The 
sum  of  these  terms  we  will  call  7 /.  As  the  physical  size  of  the  particle  changes  from  small 
to  large,  7  must  go  smoothly  from  some  large  value  to  7/.  This  transition  is  most  sensitive, 
and  hence  best  modelled,  where  the  Rayleigh  scattering  and  absorption  are  roughly  equal. 
This  occurs  since  the  Rayleigh  absorption  term  and  the  anomalous  diffraction  absorption 
term  are  close  to  or  may  be  equal  while  the  scattering  terms  are  usually  quite  different. 
This  can  play  significantly  on  the  balance  between  ^smau  and  ^iaTge  as  the  physical  size 
varies.  This  has  been  taken  into  account  by  an  additional  term  dependent  on  the  physical 
size.  This  size  is  represented  by  the  spherical  volume  equivalent  radius  of  the  spheroid. 
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From  an  empirical  lit  on  the  above  terms,  7  becomes: 


,  (54  -  7/) 

7  “  71  +  1  +  (or»/3  6)4’ 


5  3  kc 

a  =  «  + 


kc  = 


(n2  -  l)2 


7/ 


2  ifcc  +  k’  "c  66n 
2  16u2 


,  n-1  f  4n  I2  16u2  k  .  . 

~  I+((»-lW*+ip'M[N/l6  +  n’  J  +^T7T’  **  ~  (2(n  —  1)P 


We  have  verified  that  none  of  the  above  modelling  of  7  changes  the  empirical  fact  that  B 
goes  confluent  in  the  best  fit  and  hence  [42]  still  remains  valid.  It  should  also  be  pointed 
out  that  the  bridging  function  is  not  necessary  if  C7ert  is  only  required  outside  the  transition 
region.  In  this  case  ^ small  or  Qext  =  Q larger  depending  on  the  region  of  interest. 


3.0  RESULTS 

The  complete  formula,  as  presented  in  the  previous  chapter,  gives  correct  asymptotic 
behaviour  for  both  large  and  small  |m  — 1|  and  6.  Thus,  in  studying  the  error  behavior  of  the 
approximation,  the  mid-ranges  of  jm  -  lj  and  b  are  of  greatest  interest.  In  this  chapter  we 
will  compare  the  analytic  approximation  with  the  exact  T-matrix  method  as  implemented 
by  Barber  (Ref.  6)  or,  for  efficiency,  the  Mie  theory  when  r  =  1. 

Figures  1  and  2  show  the  comparison  of  Qext  versus  b  for  aspects  2  and  1/2,  re- 
spectively.  The  refractive  index  is  n  =  1.3,  close  to  that  of  water.  It  is  clear  that  the  error 
decreases  at  either  extreme  of  b.  (The  deviation  seen  in  Fig.  1  for  b  >  23  is  caused  by  ill 
|  conditioning  in  the  T-matrix  code).  The  largest  errors  are  near  and  around  the  first  two 

peaks.  This  occurs  since  much  more  of  the  scattering  physics  must  be  considered. 
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FIGURE  1  -  Comparison  between  approximation  and  T-matrix  method  for  an  in¬ 
dex  of  1.3  and  an  aspect  ratio  of  2 


FIGURE  2  -  Comparison  between  approximation  and  T-matrix  method  for  an  in¬ 
dex  of  1.3  and  an  aspect  ratio  of  1/2 
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Figures  3  and  4  are  contour  graphs  of  the  percent  error  between  the  analytic  ap¬ 
proximation  and  the  Tmatrix  computation.  The  refractive  index  varies  as  1  <  n  <  3  and 
10~5  <  k  <  10.  Figure  3  is  for  a  sphere  and  Fig.  4  is  for  a  prolate  spheroid  of  aspect  2. 
There  are  three  features  of  note.  One  is  the  increase  of  the  error  for  small  k  and  large  n. 
This  error  is  shown  well  in  Fig.  5  for  the  case  of  a  sphere  and  index  m  =  1.8  -  Oi.  This 
feature  is  simply  due  to  resonant  surface  waves  that  are  not  modelled.  Note  that  for  aspect 
2,  these  errors  are  significantly  smaller  since  the  surface  waves  are  damped  by  the  aspheric- 
ity.  The  remaining  errors  arise  from  inaccurate  modelling  of  the  bridging  region  —  that  is 
between  the  Rayleigh  region  and  the  first  peak.  The  second  feature  occurs  approximately 
when  2  <  k  <  10.  Here,  for  large  particles,  the  coherence  effect  of  the  internal  refracted 
wave  is  significant  and  hence  has  not  been  properly  modelled  (Ref.  15).  This  is  shown  in 
Fig.  6  for  an  oblate  particle  with  an  aspect  ratio  of  1/2  and  m  =  1  —  3i.  This  becomes 
insignificant  for  larger  k  since  the  particle  becomes  reflective.  The  remaining  error  at  these 
large  values  of  k  and  6  «  1  is  again  due  to  difficulties  in  the  bridging  function  attempting 
to  model  the  electromagnetic  field  on  the  surface.  The  third  significant  feature  occurs  when 
the  Rayleigh  scattering  and  absorption  are  roughly  equal.  This  occurs  since  the  Rayleigh 
absorption  term  and  the  anomalous  diffraction  absorption  term  are  close  to  or  may  be  equal 
while  the  scattering  terms  are  usually  quite  different.  This  was  included  in  the  modelling 
of  7  in  the  previous  chapter  by  kc,  [43].  Residual  errors  can  still  be  seen  in  Figs.  3  and  4 
due  to  imperfect  modelling.  These  errors  follow  curves  with  k  <x  (n2  —  1  )2/n. 


Log  Imaginary  Index  w  Log  Imaginary  Index 
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Max  Rel  Error,  Analytic  vs  T-Mat,  r=1 


E  3  -  Maximum  relative  error,  in  percent,  between  analytic  approximation 
and  Mie  theory,  r  =  1 


Max  Rel  Error,  Analytic  vs  T-Mat,  r-2 


FIGURE  4  -  Maximum  relative  error,  in  percent,  between  analytic  approximation 
and  T-matrix  r  =  2 


Qext  Qext 
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-  Comparison  between  approximation  and  Mie  theory  for  an  index  of 
1.8  and  an  aspect  ratio  of  1;  significant  surface  waves. 


FIGURE  6  -  Comparison  between  approximation  and  T-matrix  method  for  an  in¬ 
dex  of  1  —  3i  and  am  aspect  ratio  of  1/2;  significant  interna]  wave 
coherence  effect. 


Qext  _  Qext 
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Semi-Axis  b 


-  Comparison  between  approximation  and  T-matrix  method  for  an  in¬ 
dex  of  3  -  Oi  and  an  aspect  ratio  of  2;  incipient  MDR  at  b  ss  1 


Semi-Axis  b 

FIGURE  8  -  As  Fig.  7  but  T-matrix  results  low  pass  filtered  for  b  >  1 
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For  large  n,  body  resonances  can  occur.  These  are  sometimes  called  moiphology 
dependent  resonances  (MDR).  For  spheres  these  occur  near  nx  =  hr,  where  l  is  a  natural 
number.  Figure  7  shows  an  example  of  an  incipient  MDR  on  the  first  diffraction  peak  of 
Qext-  Here,  m  =  3  —  Oi  and  r  =  2.  Note  that,  despite  the  significant  perturbation  in 
the  transition  region  due  to  the  MDR,  the  approximation  is  excellent.  Only  the  first,  and 
hence  simplest  MDRs  are  modelled  (by  z\,  [5]  and  7?,  [6]).  To  show  the  accuracy  of  the 
approximation,  the  T-matrix  was  low  pass  filtered  for  6  >  1  and  is  graphed  in  Fig.  8.  Since 
by  far  the  main  contribution  is  now  from  the  diffraction  peaks,  the  underlying  accuracy  of 
the  approximation  is  apparent. 

The  next  example  is  a  model  of  extinction  by  randomly  oriented  copper  flakes  in 
the  infrared  (m  =  35  -  35  i).  An  oblate  spheroid  with  an  aspect  rttio  of  0.333  was  used.  For 
this  index,  lower  aspect  ratios  could  not  be  considered  since  the  T-matrix  will  not  produce 
usable  results  and  hence  no  comparison  could  be  made.  Qext  for  this  case  is  shown  in  Fig. 
9.  Caution  is  required  since  the  T-matrix  for  b  >  .25  begins  to  decrease  rapidly  and  will 
go  negative  for  larger  b.  This  shows  an  advantage  of  the  approximation.  It  can  estimate 
Qext  for  combinations  of  n,  k,  b  and  r  when  the  T-matrix  cannot.  Figures  10  and  11  are 
examples  of  Qext  for  water  prolate  and  oblates  spheroids  at  9.4  GHz,  respectively.  Both  are 
difficult  cases  since  the  T-matrix  is  almost  ill  conditioned  and  several  MDRs  are  becoming 
apparent.  In  Fig.  10,  the  first  MDR  is  reasonably  well  modelled  while  the  second  is  not. 
However,  the  latter  only  introduces  an  error  of  about  20  %.  Again  Fig.  11,  the  first  MDR 
is  well  modelled  and  the  second  is  not.  In  contrast  to  Fig.  10  it  introduces  a  smaller  error 
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FIGURE  9  -  Comparison  between  approximation  and  T-matrix  method  for  an  in¬ 
dex  of  35  -  35  i  and  an  aspect  ratio  of  1/3 

since  it  is  coincident  with  the  first  diffraction  peak.  Note  that  at  the  highest  b  shown  the 
T-matrix  has  become  ill  conditioned. 

Extensive  computations  have  been  carried  out  to  indicate  the  error  of  the  approx¬ 
imation  over  the  complete  range  of  stability  of  the  T-matrix  method.  Both  the  error  and 
location  of  the  error  as  a  function  of  refractive  index  m  and  aspect  ratio  r  are  given  in  the 
Appendix. 

All  our  approximate  Q ^  diagrams  in  this  chapter  and  the  appendix  were  produced 
at  a  rate  of  greater  than  104  times  faster  than  by  the  T-matrix  code.  Since  the  T-matrix 
scales  as  at  least  the  cube  of  the  optical  size  whereas  the  analytic  approximation  is  optical 
size  independent,  larger  size  parameters  or  larger  refractive  indices  lead  to  larger  speed-up 


Qext 
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FIGURE  10  -  Comparison  between  approximation  and  T-matrix  method  for  an  in 
dex  of  8.075  -  1.824  i  and  an  aspect  ratio  of  2;  water  at  9.4  GHz. 


FIGURE  11  -  Comparison  between  approximation  and  T-matrix  method  for  an  in¬ 
dex  of  8.075  -  1.824  i  and  an  aspect  ratio  of  1/2;  water  at  9.4  GHz. 
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factors.  For  example,  letting  m  =  1.5  and  r  =  2,  the  following  execution  times  on  an  Intel 
i860  40  MHz  coprocessor  are  found:  0.25  s,  b  =  1,  4.07  s,  b  =  10  and  2290  s,  b  =  75.  Other 
values  of  m  or  r  are  likely  to  increase  the  time.  The  computation  time  of  the  approximation 
is  nearly  independent  of  optical  size  and  aspect  ratio  and,  for  the  same  coprocessor,  the 
execution  time  is  .0029  s.  The  savings  in  time  is  obvious  for  large  b  and  is  still  significant 
for  smaller  b  if  many  values  of  Qext  are  required. 

4.0  CONCLUSIONS  AND  LIMITATIONS 

We  have  presented  an  approximation  to  Qext  for  randomly  oriented  spheroids  which 
gives  good  results  with  little  loss  in  accuracy  for  all  size  parameters,  aspect  ratios  and 
refractive  indices  n  >  1  and  k  >  0  and  /x  =  1.  If  high  precision  is  not  required,  the  formula 
is  far  more  economical  in  computer  time  than  the  T-matrix  method  for  obtaining  Qext- 
When  both  the  range  of  demonstrated  validity  and  the  accuracy  are  taken  into  account, 
this  formula  is  superior  to  all  other  approximations  known  by  the  authors. 

Extensive  computations  have  been  carried  out  to  indicate  the  error  of  the  approx¬ 
imation  over  the  complete  range  of  stability  of  the  T-matrix  method.  Both  the  error  and 
location  of  the  error  as  a  function  of  refractive  index  m  and  aspect  ratio  r  are  given  in  the 
appendix. 

Several  limitations  of  the  previous  numerical  approach  (Ref.  1)  have  been  removed. 
The  remaining  limitations,  that  for  n  <  1  and/or  k  <  0,  are  not  modelled  since  new  and 
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significant  physical  phenomena  arise  (e.g.  total  internal  reflections  and  optical  gain).  Even 
n  modestly  less  than  1  can  cause  problems.  Another  limitation  occurs  for  2  <  k  <  10 
and  large  particles.  In  this  region,  the  absorption  is  not  well  modelled.  If  this  effect  was 
properly  accounted  for,  Qab*  and  hence  Qtca  could  be  globally  and  readily  obtained  by  using 
the  same  approach. 


The  authors  would  like  to  thank  Dr.  Peter  Barber  of  the  Clarkson  College  of  Tech¬ 
nology  for  generously  making  available  his  T-matrix  spheroid  code  without  which  this  work 
could  not  have  been  carried  out. 
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APPENDIX  A 

Survey  of  Maximum  Relative  Error  and  its  Location 

An  extensive  computation  of  Qat  for  1  <  n  <  2,  10-5  <  k  <  10,  1/5  <  r  <  5  and 
1/10  <  rb  <  30  was  performed  using  the  T-matrix  and  the  analytic  approximation.  When 
r  =  1  a  Mie  code  was  used  instead.  From  this  calculation  the  maximum  relative  error  can 
be  deduced.  In  cases  were  the  error  may  occur  outside  the  above  limits  (such  as  values  of  m 
near  1),  but  still  inside  the  stable  region  of  the  T-matrix,  the  calculation  was  extended  to 
cover  this  region.  In  addition,  the  value  of  b  were  this  maximum  occured  was  also  recorded. 

The  sequence  of  following  diagrams  gives  the  results  of  this  computation.  The 
upper  diagram  is  the  maximum  relative  error  in  percent  and  the  lower  diagram  is  the 
location  of  the  maximum.  Care  was  taken  to  insure  that  the  T-matrix  was  stable.  This  is 
done  by  verifying  the  convergence  of  the  T-matrix  itself  where  possible,  by  comparing  the 
backscatter  efficiency  with  the  Fresnel  coefficient  estimate,  and  with  the  Rayleigh  formula. 
Since  it  is  difficult  estimate  the  accuracy  near  the  limit  of  stability  of  the  T-matrix  code 


some  additional  error  due  to  the  T-matrix  code  could  be  present. 
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